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Abstract 

We consider the precise quantum state of two trapped, coupled Bose Ein- 
stein condensates in the two- mode approximation. We seek a representation 
of the state in terms of a Wigner-like distribution on the two-mode Bloch 
sphere. The problem is solved using a self-consistent rotation of the unknown 
state to the south pole of the sphere. The two-mode Hamiltonian is pro- 
jected onto the harmonic oscillator phase plane, where it can be solved by 
standard techniques. Our results show how the number of atoms in each 
trap and the squeezing in the number difference depend on the physical pa- 
rameters. Considering negative scattering lengths, we show that there is a 
regime of squeezing in the relative phase of the condensates which occurs for 
weaker interactions than the superposition states found by Cirac et al ( quant 



ph/9706034| , 13 June 1997). The phase squeezing is also apparent in mildly 



asymmetric trap configurations. 
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I. INTRODUCTION 



Traditionally, a Bose-Einstein condensate (BEC) is often viewed as a coherent state of the 
atomic field with a definite phase . It is well known that there are problems with this view, 
however, related to the fact that the phase of the atomic field is not an observable |l|-f§]. 
The Hamiltonian for the atomic field is independent of the condensate phase and so the 
correct coherent state is only defined up to its mean number. Often it is convenient to 
invoke a symmetry breaking Bogoliubov field to select a particular phase, but this does 
not correspond to any physical field so the procedure is not totally satisfactory in a formal 
sense. In addition, a coherent state implies a superposition of number states, whereas in 
the current single trap experiments there is a fixed number of atoms in the trap (even 
if we are ignorant of that number), and the state of a single trapped condensate must be 
a number state (or more precisely, a mixture of number states). Both these problems are 
bypassed by considering a system of two condensates for which the total number of atoms 
N is fixed. Then, a general state of the system is a superposition of number difference states 
of the form 

N 

| ) = Y,c k \k,N-k). (1) 

fc=0 

As we now have a well-defined superposition state, we can legitimately consider the relative 
phase of the two condensates, which is an Hermitian observable. Indeed, the dramatic 
observation of interference between two coherent BEC's || constitutes a measurement of 
exactly this. In the absence of atomic collisions, the expansion coefficients in Eq. (H) obey 
a binomial rather than Poissonian distribution as would be expected for a coherent state. 

However, there is a more straightforward objection to the identification of the conden- 
sate with a coherent state. This is that in real experiments, the atoms experience collisions 
introducing a nonlinearity into the Hamiltonian for the system. We then know immediately 
that (unless very strongly damped) the true state can not be a coherent state, leaving aside 
issues of absolute versus relative phase. In the first treatment of this question, Lewenstein 
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and You JT0| suggested the condensate is actually in an amplitude quadrature eigenstate. In 



a fuller analysis, Dunningham et al flT| , |T2|| have shown that for positive (repulsive) interac- 
tions the state is strongly number squeezed, and resembles a bent version of the amplitude 
squeezed state that minimizes number fluctuations. This has the potentially observable con- 
sequence of increasing the revival time in collapses and revivals of the relative phase [OHTB 



due to the reduced number variance of the squeezed state. The approach of Dunningham 
et al is based on the symmetry breaking picture described above. Their model thus de- 
scribes the quantum state of a single damped driven condensate with the phase determined 
by some much larger reference condensate which does not appear in the calculation. Thus 
while the number squeezing they predict is intuitively natural, the model faces the same 
formal difficulties mentioned above in relation to symmetry breaking. 

In this paper, we combine these two ideas by seeking an accurate description of the 
ground state beyond the coherent state picture, for a system of two coupled condensates 
with a fixed total number of atoms. We do this by reducing the full quantum field theoretical 
description to an approximate two-mode problem, valid for condensates of a few thousand 
atoms. The problem is then well-defined in the senses discussed above: we deal with relative 
rather than absolute phases, and are able to consider a completely closed system without 
the complications of driving and damping, so that the ground state is unambiguously a pure 
state. Using a variational approach, we then find approximate solutions to the two-mode 
problem which are natural analogues of the single condensate states found by Dunningham 
et al. Our approach also works for negative (attractive) interactions and we predict a 
regime of significant phase squeezing in between the coherent state-like behavior with no 
interactions, and the Schrodinger cat states reported previously |i7|JT8|] that occur with 
significant interactions. 

The paper is structured as follows. In section [TI| we briefly summarize the quantum field 
theory for the two-condensate problem and derive the approximate two-mode Hamiltonian. 
In section |JTI| we discuss a representation of the two mode states using the Bloch sphere and 
outline our method for finding the ground state of the system. We construct the solution in 



detail in section [TV|. In section [V], we present our results and compare the predictions of our 
method with exact solutions for systems with small numbers of atoms. We consider negative 
scattering lengths and the associated phase squeezing in section [VT| before we conclude. 



II. THE GOVERNING HAMILTONIAN 



A. Reduction to two-mode Hamiltonian 



Our model describes two condensates of atoms of mass m, with a linear Josephson 
coupling and weak nonlinear interactions. We consider a single trap with the condensates 



distinguished by their internal atomic state fll? ]. The coupling is provided by laser-induced 
Raman transitions between the two atomic states. Following Cirac et al [pLT|| , the second 
quantized Hamiltonian takes the form 



H = H\ + H2 + H u 



H 



coup; 



with 



d 3 x^j(x) 



■ 5S v> + V5(x) + - sr ^J (!t) ^ M 



ti ( x ) 



H r 



d 3 x^J (x) ip\ (x) (x) ^ 2 (x) 



H 
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2m 

- I d 3 x 



i>x (x)^ (x)e- i5t + Vl (x) V>2 (x)e 



(2) 

(3) 

(4) 
(5) 



where j = 1, 2. Here, the field operators ipi (x) and ip2 ( x ) annihilate atoms at position x 



in condensates 1 and 2 respectively, and satisfy the relation ^ (x) , (x' 



5ij5 (x — x' 



The term Hi 2 describes each of the condensates in the absence of interactions with the 
other. They experience spherical harmonic trap potentials V\ t 2 of frequency uj\ and co>2, 
and have scattering lengths ai and a>2 respectively. The cross-phase modulation term H int 
describes collisional interactions between the condensates with scattering length a 12 . The 
laser induced coupling is described by H coup with Q the Rabi frequency and 5 the detuning 
of the classical laser field. In our work, we assume equal scattering lengths a x = 02 = a^, 
but allow the trap frequencies to differ. 



The procedure to obtain the two- mode Hamiltonian is well-known fll7|J19|| . We approxi- 
mate the field operators as ipi (x) = b\<pi (x) and ip 2 ( x ) = ^2^2 ( x ) where 4>i i2 (x) are (real) 
normalized mode functions for the two condensates, and b\ 2 are the associated mode anni- 



hilation operators which obey the standard commutation relations [bi,bj] = 0, 
Then Eq. (□) becomes 



b- fct 



Si 



Hfx(w 1 + 5) b% + u 2 b\b 2 + xAblhh + X2b\b\b 2 b 2 + Xnb\b x b 2 b 2 - ^ (b x b\ + b\b 2 ) , (6) 



where 



UJi 

Xi 
X12 



2 

U12 
2 

n 



d 3 r fa (r) 

d 3 r 
d 3 r 
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02 (r) 



77 = — / d 3 r0i (r)0 2 (r) 



(7) 

(8) 
(9) 
(10) 



for z = 1,2. Here we have moved into the interaction picture and introduced dimensionless 
variables, scaling the Hamiltonian by an appropriate energy Hujq and the position by the 
length scale Xo = Jhj muj® such that r = x/xo and (pi (r) = x^ 2 (f> (x). Further, A, = uJi/ujQ 
and Ui = 4nai/xo, U\ 2 = 47iai 2 /x . 

As shown by Milburn et al JOJ , the same two mode model also describes coupling between 
condensates in a double well potential. In this case, the two lowest modes are strictly 
speaking the symmetric and antisymmetric modes of the entire double well, but for a large 
dividing potential barrier it is an accurate approximation to use modes describing atoms in 
one or the other trap. The linear coupling is provided directly by spatial tunnelling through 
the barrier and has a strength 77 = AE/ujq where AE is the frequency separation of the two 
(linear) modes |19|]. Assuming the potential barrier is relatively strong, the modes are well 
separated and we can neglect H int . 

Equation defines our problem completely. For large condensates, the mode functions 
are altered by the collisional interactions and the two-mode approximation breaks down. 
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As shown in Ref. ||19|| , a simple estimate shows this occurs when the number of atoms N, 
satisfies Na ^> x , where a is a typical scattering length and x is a measure of the trap 
size. Assuming a « 5 nm |],[2(J and a large trap with xq ~ 10 fim, we find the two-mode 
approximation should be acceptable for iV < 2000. 

B. Angular momentum representation 

The Hamiltonian can be reduced to a simpler form by exploiting the equivalence between 
the algebra for two harmonic oscillators and that for angular momentum, by introducing 



the new operators 21,19 



J+ = b\b 2 , 



J- = hbl (11) 
J z = \ (b\h - b% 



and 

J.= \ (J+ + J-) , (12) 
J y =^-.(J + -J-), (13) 

These operators do indeed satisfy the usual angular momentum commutation relations, 
justifying the choice of notation. In addition we find 

J 2 = Jl + Jy +Jz = y{j + 

where the total number of atoms N = n\ + n 2 = b\bi + 6^2 is a constant of the motion, so 
we deduce that we are working with the angular momentum algebra for J = N/2. In terms 
of the new variables, the Eq. (|6|) takes the form 

H = J (u>i + 6 + lj 2 - xi - Xi) + J 2 (Xi + X2 + Xn) + AljJ z + X+ J 2 Z - | ( J+ + J-) 
= AljJ z + x+Jl - vJx, (15) 
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where in the last line we have dropped an unimportant constant, and we have introduced 
the effective detuning 

Alj = lu 1 + 5-u 2 + (2J-1)x-, (16) 

and effective nonlinearity 

X+ = Xi + X2 - Xu- (17) 

Note the useful fact that both the difference in the self-nonlinearities X- = Xi ~ X2 and 
the cross-nonlinearity X12 merely shift the values of Au and x+ an d introduce no new 
terms. Thus there is no restriction of the physics by assuming equal scattering lengths. We 
now calculate these parameters for realistic experimental values. Taking the Na 23 atom for 
example, we have a pa 5 nm, and suppose trap frequencies of order Ui = 1000 s _1 ||. Then 
taking the scaling frequency uq = 1 s _1 , we obtain uoi pa 1500, and x% ~ 1-4- Therefore, the 
detuning Au may range from zero to a few hundred. The coupling strength rj is largely 
arbitrary. In the spatial case, it can take values up to the order of the trap frequencies |T9| , 
or can be made as small as desired by incresing the trap separation. 

III. OUTLINE OF APPROACH 

For the remainder of the paper we are concerned with the ground state of Eq. (|15|). 



Milburn et al [IS] have presented numerical calculations of the energy spectrum of this 



Hamiltonian and the dynamical problem has also been studied |19| , |22| . Rather than the 
spectrum or dynamics, however, our concern is with the detailed properties of the lowest 
eigenstate and their dependence on the effective detuning and nonlinearity. In general, the 
eigenstates of Eq. ( |i5|) can not be written analytically. For systems with at most a few 
hundred atoms, it is feasible to find the exact eigenstates in the basis of J z eigenstates 
I J, m = — J, . . . , J) z numerically. Our semi-analytic approach can be used for systems of 
arbitrary size and lends considerable insight to the problem. 
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A. Bloch sphere 



Our approach relies closely on the Bloch sphere representation of angular momentum 



which we must briefly introduce. A detailed analysis has been given by Arecchi et al [23 
Quantum states in the angular momentum Hilbert space can be usefully represented on the 
Bloch sphere. Certain states — the atomic coherent states or "Bloch" states [^31 — correspond 
to a single point on the sphere. Defined as the rotated states \9, tp) = Re tV \J, —J) z , where 
the rotation operator 

Re ;V = exp [0/2 (j+e~^ - J_e^)] = exp [-id {J x sin <p - J y cos tp)] , (18) 

they are labeled by the spherical coordinates 9 and p corresponding to the state's point 
on the sphere. Note that in terms of our BEC problem, the north pole \9 = n) and south 
pole \9 = 0) represent the states with all atoms in mode 1 or 2 respectively. States lying on 
the equator with 9 = tt/2 represent an equal division of atoms between the modes, (which 
for 77 7^ 0, does not imply the number state \N/2, N/2) , but rather an entanglement of the 
form (JJ) with a binomial distribution of expansion coefficients). The Bloch states are the 
analogs in the angular momentum algebra of the standard coherent states of the harmonic 
oscillator [J£J. They share a number of properties with the coherent states, for instance 
minimum uncertainty in the natural variables. In addition, more general non-classical states 
described by the state vector \ip) or density matrix p can be naturally pictured in terms of 
a quasi-probability distribution function on the sphere 

Q m (9,p) = \(9,p\iP)\ 2 , (19) 
Q p (9,p) =(9,p\p\9,p), (20) 



in analogy to the familiar Q function in the harmonic oscillator phase-plane P3 |. Func- 



tions analogous to the standard Glauber-Sudarshan P and Wigner distributions can also be 
defined. As discussed below, these analogies can be made precise using a formal contrac- 
tion from the angular momentum Hilbert space to the Hilbert space for a single harmonic 
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oscillator |23|. While the Bloch states lack some of the useful properties of the coherent 



states [^,^5|], nonetheless, we see below that the Wigner or Q functions make for useful 
measures of quantities such as the squeezing in the number difference or relative phase in 
the ground state. 

B. Mathematical procedure 

The angular momentum commutation relations make a direct solution to our problem 
in the full Hilbert space difficult. The previous section suggests the following alternative 
approach. We assume the ground state we seek has a quasi-probability distribution localized 
to a particular part of the Bloch sphere. (Thus we immediately exclude Schrodinger cat 



states such as those found by Cirac et al |17| and Ruostekoski et al ||18|| , which we treat 
numerically in section |VT| .) We apply a rotation to the Hamiltonian to bring the mean value 
of the state to the south pole of the Bloch sphere. This must be done self-consistently as we 
do not actually know the mean value of the state until we have solved the problem. We then 
project the problem to the harmonic oscillator phase plane using the contraction operation 
of Ref . [^| . The problem can then be solved in the plane and the ground state plotted as a 
Wigner distribution. Finally, we project this distribution back on to the sphere, and rotate 
it to the original mean value which has by now been determined. Equivalently, we can think 
of the problem being solved in the oscillator phase plane that is tangent to the sphere at the 



mean value of the state. In section IV we make these ideas precise and provide the solution. 



C. Linear problem 

In the absence of the nonlinear term (x+ = 0), the problem is trivial. The Hamilto- 
nian becomes 

H = AljJ z - r}J x , (21) 
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and using the rotation relations in appendix |A], it is easy to see that the rotated Hamiltonian 
H' = Rq^HRqI. where tan6> = 77/ Au>, is just a multiple of J z with ground state \J, —J) z ■ 
Inverting the rotation, the exact ground state for the original Hamiltonian is simply the 
Bloch state \9,<p) = |tan _1 (rj/Au) , 0) . Thus the the ground state of coupled ideal gas 
condensates is the entangled state analog of the coherent state. This result is well-known, 



though it is more commonly expressed in the number difference basis ||26|| . As expected, 
for vanishing Au, the traps are equivalent and there is an equal number of atoms in each, 
whereas for AlD 7^ 0, the ground state has more atoms in the weaker trap. Note that the 
relative phase of the condensates 7? = 0. This is the reason for our choice of the negative sign 
in front of 77 in the original Hamiltonian (pl|) . Even when we consider the more complicated 
states of the nonlinear system, we see by symmetry that the mean value of the state must 
still have ip = 0, simplifying the rotation operators we need to consider. In practice, the 
phase of 77 is determined by the phase of the driving laser field. By a suitable choice of 
coordinates we may always take it to be zero. 



IV. NONLINEAR PROBLEM 

A. Contraction from angular momentum to harmonic oscillator Hilbert space 

The nonlinear problem is much more involved. As indicated earlier, the first step is 
to perform a rotation of the Hamiltonian by an undertermined angle 9 and then project 
the new Hamiltonian into the harmonic oscillator phase plane with operators a, satisfying 
1. This procedure can be made rigorous through the concept of a group contraction 
from the angular momentum Hilbert space to the harmonic oscillator Hilbert space. The 
details can be found in Ref. p3[ to which we refer the interested reader. Quoting the results, 



a, 



the contraction is made by the identification of operators according to 



J + - l -a\ (22) 
J_ -> -a, (23) 
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J z -> tfa - — , (24) 

where c = l/\/2J and the spaces are formally identical in the limit c — ► 0. In the same 
limit, we can contract eigenstates of J z to the harmonic oscillator number states: 

\J,M) -> |n = J + M), (25) 

and we relate the coordinates according to 

- exp (i<£>) — > ca. (26) 

Later we also use the quadrature operators X = a + a) and Y — —i (a — aM , which are 
the contractions of J x and respectively. Geometrically, we visualize this contraction as 
a projection from the Bloch sphere to the phase plane with the south pole of the sphere 
coincident with the origin of the phase plane. Note that the coherent amplitude a can take 
values throughout the whole phase plane only in the limit c — > and there is naturally a 
distortion involved in the projection. Physically, by performing the contraction we discard 
the knowledge that the true ladder of states is bounded at both ends rather than just the 
lower end. However, providing the state is localized near the south pole and c < 1 (large 
atom number), the distortion is small. The contraction process also maps functions from the 
sphere to the plane, so for instance we can identify the (rotated) Bloch sphere distribution 
function Q (6, ip) with the standard phase plane function Q (a) = (a \p\ a). Here, we a define 
the Wigner-like distribution on the sphere by a projection of the harmonic oscillator Wigner 
distribution using Eq. (|26|) in reverse. 

B. Variational solution 

1. Gaussian part 

We are at last ready to find our approximate solution to the full problem. The procedure 
is somewhat involved mathematically, and we state only the main intermediate steps here, 
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leaving the details to appendix The first step is to rotate the Hamiltonian (|I5D around 
the positive y-axis by an undetermined angle 9, and perform the contraction operation to 
the single oscillator phase space to find the new Hamiltonian F. Following Dunningham et 



al ||12||, we write this Hamiltonian as 



Fg + Fng, 



(27) 



in terms of a Gaussian part Fg and non-Gaussian part F^g, the latter of which satisfies the 
constraints 



(F NG ) = 0, 
([a,F NG ]) = ([a\F NG 
([a, [a,F NG ]]) = (\a, \a\F NG 



(28) 



0. 



This separation allows us to find the ground state of the Gaussian part first, and by assuming 
a weak nonlinearity, treat the non-Gaussian part as a perturbation. Appendix [A] contains 
the expression for the non-Gaussian part. The Gaussian part is 



F G = K + L (a + a f ) + Sa ] a + T (a 2 + 



,12 



(29) 



where 



K = x+J/2 sin 2 9 + x+J 2 cos 2 9 - J (Alj cos9 + t] sin 9) 

+2v / 2Jx+ cos 9 sin^ (aV) + X + cos 2 9 ( (a) 2 a 2 ) - 2 (a 2 ) 2 - 4 (a f a 

L = ^ { Aid sin 6» - ?7 cos 6» + x+ sin 9 cos [- (2 J - 1) + 2 <V J) + 4 (Va)] } 
S = r? sin 6» + Au cos 6> + x+ {j sin 2 9 + cos 2 9 [- (2 J - 1) + 4 (a f a)] } , 



T = x+ ( ^ sin 2 + cos 2 (a 2 



(30) 

(31) 
(32) 
(33) 

In Eqs. (p0|)-(p3|), we have taken (a 2 ) = (<^ 2 ) which follows from the choice of rj as real. 
Note that Fq depends on moments taken over the state which is the solution we are seeking. 
We can solve the Gaussian part to different levels of accuracy according to how we account 
for these expectation values. 



12 



a. Non self- consistent approach We first assume we have known values for the expecta- 
tion values. For example, we may take a mean-field approximation in which all the moments 
are zero, or as explained below we may have obtained estimates for the moments from a 
previous less accurate calculation (such as the mean- field one). We consider a self-consistent 
approach in the next section. The rotation angle 9 is fixed by requiring that the linear terms 
should vanish: 

L = 0. (34) 

Except for a constant term, Fq is now purely quadratic and we perform a Boguliobov 
diagonalization by writing 

a = b cosh r — sinh r, 

a) = cosh r — b sinh r. (35) 
Substituting Eqs. (^) into Eq. (|29| ) and setting the terms in b 2 and tf 2 to zero, we obtain 



/ Q — OT 

exp( - 2r) = VsT2T' (36) 



while the diagonalized Hamiltonian is 



Fq = K + (S cosh 2 r-2T cosh r sinh r) + VS 2 - AT 2 b ] b. (37) 

The first two terms are constants, so the ground state is just the vacuum in the b rep- 
resentation. As the transformation (|35|) is induced by the squeezing operator S (r) = 
exp r (a 2 — a) 2 ^j /2 [24]], the b eigenstates |z) transform back as 

~jj =S{r)\i). (38) 

Thus the ground state in the a representation is just the squeezed vacuum with 



X 2 ) = exp (_2r) = J|rrJp (39) 



In the mean-field limit we have the simple results 
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Alj = r] cot 9 + x+ cos 8 (2J — 1) , 



(40) 




rj + Nx+ sin 3 9 



(41) 



For Au = 0, we find symmetric states with 9 = it/ 2 as is natural. Moreover, in the limit 
9 — > 7r/2 (that is, for Au/ [(2 J — 1) x+] — ► 0), the non-Gaussian part of the Hamiltonian 
Fng = (see appendix [A]) and Fq is independent of any expectation values. Thus in this 



limit, the projected state is exactly a squeezed state with (X 2 ) = yr// (rj + x + iV).We note 

in passing that for a negative nonlinearity, Eq. ([|l]) predicts (X 2 ) > 1 which indicates a 
possibility of phase-squeezing. We return to this in section [V]]. 

b. Self- consistent Approach We can also find the ground state of Eq. (|29|) with a self- 
consistent approach in which the expectation values are determined to Gaussian approxi- 
mation in the course of the calculation. In this case, the first two terms in Eq. (^) can 
not be considered constants and in general the 6-vacuum is not the lowest eigenstate. The 
correct approach is to assume a squeezed vacuum solution |r) = S (r) |0) to Eq. ( |29| ) and 
find 9 and r by minimizing the expectation value of the energy subject to the constraint 
in Eq. (0). Performing the transformation (|35|) with this value of r, gives a Hamiltonian 
in the b representation with a small off-diagonal part which can be transfered to the non- 
Gaussian part Fng yet to be treated. In fact, we have found that we can obtain virtually 
identical results by proceeding directly with the Boguliobov diagonalization, and solving 
Eqs. (|34] ) and ([39]) simultaneously for 9 and r, where L, S and T now depend on r through 
the quadratic moments. 



We now include the effects of the non-Gaussian part F^g as a perturbation to the 
squeezed state just found. We use second-order perturbation theory to write the corrected 
state as 




2. N on- Gaussian part 




(42) 
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where from Eq. Q B^ 1 = AVS 2 - 4T 2 is the energy of the unperturbed state 
then calculate the Wigner distribution 




W(a) 



d ze 



(44) 



where the symmetric characteristic function is x( z ) = Trjpa ex P ( zc ^ ~~ z " 



)} and p a = 

. The details of these calculations are given in appendix ||. The final answer 
has a closed form in terms of Hermite Gaussians but is too long to write here. 



3. Combined Approach 

In finding ground states we have used the above steps in an iterative scheme. A first 
approximation is found using the self-consistent approach to the Gaussian part followed by 
the perturbation theory. The quadratic moments appearing in Eqs. (0) and (|39| ) [through 
L, S and T] are calculated using this first approximation and a new Gaussian state chosen. 
Finally the perturbation theory is applied again. 



4- Anti-rotation 

To complete the problem, the contours of the Wigner function just obtained must be 
projected back to the sphere and the original rotation of the Hamiltonian by angle 9 reversed. 
This is completely elementary and we reserve the equations for appendix |C|. 



V. RESULTS 
A. Exact States 



In this section, we present a mixture of exact numerical results and those obtained by 
our semi-analytic procedure. This allows us to test the agreement in the regime where the 
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size of the Hilbert space is small enough to permit a complete numerical solution. To begin, 
in Figure [I], we show three sample exact states for N = 100 atoms plotted as contours 
of Wigner functions on the surface of the Bloch sphere. There are two contours for each 
state at heights e -1 and e _1 /4 of the maximum of the Wigner function. States (a) and 
(b) show states with a nonlinearity x+ = 0-75 and detuning Au = and 30, respectively. 
Both states show strong squeezing in the number difference (the vertical axis J z ), while for 
the asymmetric case (b) the atoms are predominantly found in the trap of lower energy. 
Note that the sense of squeezing is along parallels of latitude and not along the great circle 
through the mean field point. This is rather obvious — we expect squeezing along the number 
difference axis J z , but it has the effect that the states are in most cases far from minimum 
uncertainty in the natural variables. We discuss this shortly. For comparison, the state (c) 
is just the Bloch state solution to the linear problem (x+ = 0)> with Au; = —0.44, for which 
the contours are circles. Note that the mean angle 9 is the same distance from the equator 
9 = 7r/2 for states (b) and (c) despite very different values of the magnitude of the detuning 
|Ao;| . As indicated by Eq. (]4"0"|), the positive nonlinearity tends to push states back towards 
the equator and is balanced by a much larger value of the detuning. This derives from an 
energy competition between the terms in (J z ) and (Jf) in the Hamiltonian flT5p. 



B. Comparison with the model 

We illustrate the results of our semi-analytic method by rotating the nonlinear states (a) 
and (b) in Fig. [I] to the south pole, and projecting them to the plane in Figs. ^|(a) and ^|(b) 
respectively. The contours are at e _1 /2 of the maximum of the Wigner function. Here we see 
the effect of the fact that the orientation of the squeezing is along the parallels of latitude. 
The state originally at 9 = tc/2 [Fig[l|(a)] is a precise squeezed state with no bending, but the 
asymmetric state in Fig. |l|(b) is distorted on projection [solid line in Fig. 0(b)]. The other 
lines in ^(b) indicate our semi-analytic prediction to second order perturbation theory. The 
dashed line shows the solution using the mean-field approximation for the expectation values 
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in Eqs. (pQ[) (p3|) . The dot-dashed line is for the improved result in which the expectation 
values are first estimated using the self-consistent approach. The bending we find here is 
a clear analog of that found for a single condensate by Dunningham et al |12| but in our 



case arises purely from the geometric effect of projection. As our theory gives the exact 
symmetric state, the lines are coincident in Fig.^(a). 

The dependence of the mean angular position of the state 9 = tan -1 (— ( J x ) / (J z )) on the 
detuning and nonlinearity for exact solutions with N = 200 is shown in Fig. |3](a). This, of 
course is a measure of the imbalance in the populations of each trap: (rii) = J (1 — cos#) , 
(n-z) — J (1 + cos#) . We plot the mean angle 6 as a function of the nonlinearity x+ f° r 
detunings of Alj = 0, 5, 25, 50, 75 and 100 which label the curves. As x+ increases, the 
mean value increases from the linear result rj = tan -1 (r//Au>) towards the symmetric value 
= 7r/2, with the curves for larger detuning shifting at larger nonlinearities. From Eq. (|40"|). 
we see that the most rapid change occurs for \+ ~ Au)/ (2J — 1). As explained above, the 
tendency toward symmetric states is a result of an increasing energy penalty for asymmetric 
states from the (J|) term in the Hamiltonian. We check the accuracy of our model in 
Fig. |3](b) showing the discrepancy in the mean angle 9 according to the mean-field (dotted 
line) and self-consistent predictions (solid line), from the exact value calculated numerically. 
The curves are labeled with the value of the detuning Alj. There is a clear improvement 
with the self-consistent case, though it is less dramatic for the larger detuning. 



We consider the behaviour of the spread in number difference 5n = yVar {rii — n 2 ) = 
4y / Var ( J z ) for the same parameters in Fig. |j. The solid lines are the exact result, the dotted 
lines our approximate result using the corrected quadratic moments, and again the curves 
are labeled by the detuning Au. For Alj = 0, the state is always centered on the equator 
and the number squeezing grows stronger with the nonlinearity. For this case, in the limit 
c —y when the projection gives the exact solution, we have 5n = y/N [77/ (77 + x+N)} 1 ^ ■ 
The discrepancy of this curve from the exact result is not visible in Fig. |4j. The behavior 
is somewhat different for the other cases. Initially the spread in number increases, before 
turning around and becoming coincident with the decreasing symmetric case. The initial 
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rise in the variance agrees closely with the Bloch state result 5n = \J N sin 9 (not shown in 
figure) until just before the maxima of the curves. Thus we see that initially the nonlinearity 
shifts the mean value of the state without affecting its shape. In this plot, we see that our 
approximate method is less successful for cases with large detunings. These states are highly 
asymmetric and therefore the projected states show significant bending. The perturbation 
from the Gaussian squeezed state is thus larger and our calculation less accurate. 



VI. NEGATIVE NONLINEARITIES 

We demonstrate briefly here that for a negative nonlinearity there is a regime of phase 
squeezing rather than number squeezing. Using a mean-field picture, Cirac et al [|17| have 
found a range of superposition states for negative nonlinearities (attractive interactions). 
They show the two lowest energy states are even and odd superpositions of states in which 
most of the atoms are in trap 1 or most are in trap 2. In our notation they arise as follows. 
In the mean field approximation of Eq. (f40|) and taking the symmetric case Au = 0, we have 

sin9 = ^P- -. (45) 

This equation clearly only has solutions for \x+\ sufficiently large. When this is true, there 
are two degenerate mean field ground states \9, 0) and | vr — 6,0) . Cirac et al have shown that 
the superposition or "Schrodinger" cat states |±) = ^= (\9, 0) ± \ir — 9, 0)) give a lower value 
for the energy and thus are a better approximation to the lowest energy levels. Ruostekoski 



and Walls |T8[ have proposed a scheme for generating similar superpositions in number for 
free condensates. Numerically, we have found that the exact ground states in this regime are 
indeed of a superposition nature, though of course they are superpositions of distorted Bloch 
states, not of true Bloch states. What about the regime r\ > \x+ \ (2 J — 1) for which Eq. ([45]) 
has no solutions? Equations (f40|) and ([|l|) give solutions with 9 = n/2 and (X 2 ) > 1. As 



in the Gaussian approximation the states are minimum uncertainty we are led to expect 
phase squeezing. While our method is applicable for negative nonlinearities, the states can 
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be highly non-Gaussian and the variational method is not always very successful. Therefore 
we use numerical results to indicate that the phase squeezing does indeed occur. We reduce 
the number of atoms to iV = 100 to make squeezing more obvious in the figures. Thus in 
the mean field approximation, we expect superposition states for x+ < —1/99 ~ —0.0101. 
In Figs. |5](a)-(b) we show the Wigner function for a succession of states with Aui = 0, 
and nonlinearities (a) x+ = —0.01, (b) x+ = —0.0115. For a vanishing nonlinearity we 
would have circles centered on the equator. The state becomes increasingly elongated in 
the number difference direction and strongly squeezed in the relative phase direction around 
the equator. Note that state (b) lies in the range where the mean field picture predicts a 
cat state. With further increase to x+ — —0.012, the phase squeezed state bifurcates to the 
cat state. This is seen in Fig. [|(c) where we have plotted the Q function rather than the 
Wigner function to avoid interference fringes. In Fig. |5](d) we treat an asymmetric case with 
Auj = 0.001, x+ — —0.0115. Here, the energy gained by adopting the superposition state 
is outweighed by the energy difference between the two traps and the lowest energy state is 
a single drawn out "tear-drop" . The extended tail is clearly a vestige of the superposition 
states that are favorable for vanishing or very small asymmetries. The long tail and phase 
squeezing may be thought of as a "best attempt" to attain a cat-like state. In Fig. ^|we show 
the phase variance A0 = ((^) — (Jy} 2 ^j /(J/2) as a function of the nonlinearity for several 
values of the detuning ACo. For small asymmetries there is strong phase squeezing. At larger 
detunings, the system is too far from the superposition state regime and the residual phase 
squeezing is quenched out. 



VII. CONCLUSION 

In this paper we have studied the quantum statistics of the ground state of a two-mode 
model for coupled Bose-Einstein condensates. We find strong squeezing of the number 
difference for positive nonlinearities and a regime of squeezing in the relative phase for 
negative nonlinearities. Within the validity of the two-mode approximation our model can 
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treat systems of arbitrary numbers of atoms. However, its applicability is limited by the 
eventual distortion of the condensate mode functions that occurs for condensates of more 
than a few thousand atoms. In order to treat larger condensates, one must account for a 
larger number or possibly all of the trap modes. This might be attempted by a variational 
solution of the full second-quantized Hamiltonian. In this fashion, Cirac et al [I7| have 



calculated the energies of superposition state, while Spekkens and Sipe [[FJJ have considered 



the coherence properties of double traps, but neither have discussed the detailed shape of 
the ground state. Other authors are currently using stochastic simulations of generalized 
Gross-Pitaevski equations with additional quantum noise terms to account for the higher 
modes |28[ . 
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APPENDIX A: SEPARATION OF THE CONTRACTED HAMILTONIAN 

Here we provide a fuller account of some of the steps in finding the ground state in the 
single oscillator Hilbert space. We first note that the rotation operator fll8|) transforms J x 
and J z as 

Re^JxRe^l = Jx cos 9 - J z sin 9, (Al) 
Re^JzRe^ = Jx sin 9 + J z cos 9. (A2) 

Using these relations, we rotate the original Hamiltonian ( fL5|) to obtain. 

H' = Rg t7T HRg^. 

= J x (Au) sin 9 — 7] cos 9) + J z (Aoj cos 9 + rj sin 9) 

+X+ Jx sin2 + J z cos2 6 + sin9 cos 9 (J X J Z + J Z J X ) . (A3) 
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Performing the contraction to the harmonic oscillator Hilbert space we find the new Hamil- 
tonian 

1 



F = x+J ( ^ sin 2 9 + J cos 2 9j - J(rj sin 9 + Au cos 9) 

+ (a + [Au sin ^ — 77 cos 9 — x+ si n $ cos # (2 J — 1)] 



+ (a 2 + a 1 " 2 ) xJk sin 2 + a f a {v sin + Acj cos 9 + x+ J sin 2 9 — (2J— 1) cos 2 6> } 



+ (Va 2 + at 2 a) x+v^2J sin 6 1 cos 6* + a^ 2 a 2 x + cos 2 9. 



(A4) 



Separating F into Gaussian and non-Gaussian parts by imposing the constraints in Eqs. (|28|) 
gives 

F G = x + sin 2 9 + cos 2 9 ( J 2 - (a 2 ) 2 - 2 (a f a)^ ") - J (rj sin + Alu cos 0) 

+ (a + a 1 ") ]J^{Au sin 6> - 77 cos + sin 9 cos - (2 J - 1) + 2 (a 2 ) + 4 <Ja f aJ> } 
+ (a 2 + ^ 2 ) x+ sin 2 9 + cos 2 (a 2 )) 

(A5) 



+a f a (77 sin 9 + Acu cos 6> + x+ { J sin 2 + cos 2 9 - (2 J - 1) + 4 ^a f aj> }) 



and 



-Fjvg = (Va 2 + a^aj x+v / 2J sin 6> cos 9 + a t2 a 2 % + cos 2 9 

- (a + a f ) x+ v 7 ^ sin cos 9 ((a 2 ) + 2 ^a f a^) - (a 2 + a f2 ) x+ cos 2 9 (a 2 ) 
-a ] aAx+ cos 2 9 (a) a) + x+ cos 2 9 {(a 2 ) 2 + 2 (a f a) 2 ^ . (A6) 

In these expressions, we have taken (a 2 ) = (c^ 2 ) which must be true by symmetry ((cp) = 0). 



APPENDIX B: EFFECTS OF THE NON-GAUSSIAN HAMILTONIAN 

Here we show the details of the perturbation calculation to find the effects of the non- 
Gaussian part of the Hamiltonian F^g- We show working only for the first order correction. 
The second order calculation proceeds identically but is much longer. We begin with the 
expression for the first order perturbation to the Gaussian ground state 
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It is easier to work in the b representation with the state 

K>X = |0> + Egi^|0>, (B2) 

where F^g must be expressed in the b basis. Applying the Bogoliubov transformation fl35|) 
we obtain 

F NG = x+ cos 2 9 [c 2 s 2 (V 4 + b A ) - (c 3 s + cs 3 ) (bH + tfb 3 ) + (c 4 + s A + 4cV) bH 2 ' 

+ X+ V2Jsm0cos6 [(cs 2 - c 2 s) (b ]3 + b 3 ) + (c 3 - s 3 + 2cs 2 - 2c 2 s) {b ]2 b + b ] b 2 
+ A (6 t2 + b 2 ) . 

where c = coshr and s = sinhr, and A = T (c 2 + s 2 ) — csS accounts for any quadratic part 
left over from the self-consistent approach. We have typically found this to be negligibly 
small. Substituting Eq. (|B3|) in Eq. (|B2] ) we find the unnormalized new state as 

(B4) 



(B3) 



$W) = k \0)+k 2 \2)+k 3 \3) + h\4) 



with 



k = 1, 
k 2 ~- 

h-- 

&4 - 



A 



V2VS 2 - 4T 2 

[2x+V2J sin 6 cos 6 (cs 2 — c 2 s) 
V3 V5 2 - 4T 2 

13 x+ cos 2 9c 2 s 2 
V2 v^^iP"' 



(B5) 
(B6) 

(B7) 

(B8) 



$W\ /$(!) 



Setting the density matrix p a = 

X {z) =Tr{ Pa e zat ~ z * a } 

= Tr (r) p a S (r) 5 f (r) e 2aWa S (r) } 



, we define the characteristic function 



^2 kikj (i 

{i,j}e{0,2,3,4} 



e fet(zc+2*s) e -fe(2:s+2;*c) 



■^{zc+z* s)(zs+z* c) 



(B9) 
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from which the Wigner function is found as 



W {a) = \ I e az '- a ' z X (z) tfz. (BIO) 



7T 2 



Expanding the exponential in the expectation value of Eq. (|B9| ) and using the Rodrigues' 
formula for the Hermite polynomials H n (x) = (— l) n exp (x 2 ) d n /dx n exp (— x 2 ) [2j| one finds 
that the Wigner function has a closed form expression as a sum of two-dimensional harmonic 
oscillator functions. This makes for rapid numerical calculation, but the expression is too 
lengthy to warrant inclusion. 



APPENDIX C: INVERSE ROTATION OF DISTRIBUTIONS ON THE SPHERE 

Suppose a contour Cq of the Wigner function in the plane is parametrized as 

C p (t) = (x (t),y (t)), (CI) 

where x = a + a*, y = — i (a — a*) are quadrature variables. Projecting onto the sphere 
using the inverse of Eq. ( PEP we get 

h (*) = ^xl(t)+yl(t), (C2) 

^ w - t -" 1 (Si) • < C3) 

where care must be taken in determining the correct quadrant of <^o- in Cartesian coordi- 
nates, the contour on the sphere is expressed 

Q = (j x (t) , j y (t) , j z (t)) = J (sin cos <p , sin 9 sin ip , - cos 9 ) , (C4) 

and is transformed by the rotation to 

CI = {j x cos 9 - j z sin 9, j y , j x sin 9 + j z cos 9) , (C5) 

which may then be reexpressed in terms of new spherical coordinates 9\ and Finally, if 9 
is small so that the number of atoms in trap 2 greatly exceeds that in trap 1, we can obtain 
a Wigner contour for the state of a "single" condensate by projecting the contour C\ {9±, tp\) 
directly back to the plane using Eq. (p6|). 
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FIGURES 

FIG. 1. Contours of the Wigner function on the Bloch sphere for exact solutions with N = 100 
atoms and a) detuning Ad) = 0, nonlinearity \+ = 0.75, b) Ad) = 30, x+ = 0.75, c) Ad) = —0.44, 

x+ = o. 

FIG. 2. Contours of the Wigner function projected into the plane for states (a) and (b) in 
Fig. 1. The solid lines are the exact results. In (b), we also show the prediction of the mean 
field approximation (dashed line), and that using corrected versions of the quadratic moments 
(dot-dashed line). 

FIG. 3. a)Mean angular position 9 as a function of nonlinearity \+ f° r Ad) = 0, 25, 50, 75 and 
100. b) Discrepancy in the mean angle 9 from the exact result as calculated by the mean field 
picture (dotted) and corrected moments picture (solid). The curves are labeled by the detuning 
Ad). 

FIG. 4. Spread in the number difference 5n = (Var(ni — r^)) 1 / 2 . Solid lines are exact results, 
dotted lines are predictions of the corrected moments theory. The curves are labeled by the detuning 
Ad). 

FIG. 5. Exact states for negative nonlinearities. (a) Wigner function for \+ = —0.01, Ad) = 0, 
(b) Wigner function for \+ = —0.0115, Ad) = 0, (c) Q function for x+ = —0.012, Aid = 0, (d) 
Wigner function for \+ = —0.0115, Ad) = 0.001, 

FIG. 6. Relative phase variance A<fi = (J^Jy^ — (^) 2 ) /(J/ty as a function of nonlinearity x+- 
The legend indicates line types for different detunings. 
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